Install the following packages and load them:
install.packages("devtools","expm", "ggplot2", "NPflow", "truncnorm")
devtools::install_github("borishejblum/CRPdemo")
library(NPflow)
library(CRPdemo)
library(ggplot2)
n <- 100 # datasize
set.seed(1231)
# Sample data
m <- matrix(nrow=2, ncol=4, c(-1, 1, 0, 2, 1, -2, -1, -2))
p <- c(0.2, 0.1, 0.4, 0.3) # frequence des clusters
library(expm)
s <- array(dim=c(2,2,4))
s[, ,1] <- matrix(nrow=2, ncol=2, c(0.1, 0, 0, 0.1))
s[, ,2] <- matrix(nrow=2, ncol=2, c(0.01, 0, 0, 0.1))
s[, ,3] <- matrix(nrow=2, ncol=2, c(0.1, 0.08, 0.08, 0.1))
s[, ,4] <- .1*diag(2)
c <- rep(0,n)
z <- matrix(0, nrow=2, ncol=n)
for(k in 1:n){
c[k] = which(rmultinom(n=1, size=1, prob=p)!=0)
z[,k] <- m[, c[k]] + expm::sqrtm(s[, , c[k]])%*%matrix(rnorm(2, mean = 0, sd = 1), nrow=2, ncol=1)
}
plot(t(z), pch=16, xlim=c(-3,3), ylim=c(-3,3), xlab="Dim 1", ylab="Dim 2")
We first need to define prior distribution and set hyper parameters.
hyperG0 <- list()
hyperG0[["mu"]] <- c(0,0)
hyperG0[["kappa"]] <- 1
hyperG0[["nu"]] <- 4
hyperG0[["lambda"]] <- diag(2)
# Scale parameter of DPM
alpha <- 2
#number f iteration
N <- 5 #15
gibbsDPMalgo1(z, hyperG0, alpha, niter = N, doPlot = TRUE) #, pause = 0
Try increasing the number of iterations
gibbsDPMalgo2(z, hyperG0, alpha, niter = N, doPlot = TRUE) #, pause = 0
Now let’s try to change alpha and see how it impacts the algorithm
gibbsDPMalgo2(z, hyperG0, alpha = 0.0001, niter=N, doPlot=TRUE) #, pause = 0